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I. INTRODUCTION 

The investigation of time-resolved currents in meso- 
scopic devices has gained a lot of interest over the past 
few years. This is not only because of the potential ap- 
plication to quantum computing but also due to the ad- 
vent of new experiments specifically looking into time- 
dependent electron transport P, 0. For example, ma- 
nipulation of quantum dot systems is performed by us- 
ing pump-probe schemes with a single voltage pulse. The 
rising and the falling edge of the pulse lead to pumping 
and probing the device, respectively. The experiments in- 
clude transient-current spectroscopy of single quantum- 
dots [H and coherent manipulation of charge 4] and spin 
0, @ qubits in double quantum dots (DQDs) . 

The theoretical description of the electric current 
through a device coupled to two electron reservoirs is 
usually based on Keldysh non-equilibrium Green func- 
tion (NEGF) techniques 0, [1] ■ Within this approach the 
description of piecewise constant and sinusoidal voltage 
pulses is readily possible in the wide-band limit (WBL). 
For harmonic modulations more sophisticated methods 
B im combining Floquet theory and NEGF formal- 
ism have been developed and allow going beyond WBL. 
In order to overcome the limitations of the special form 
of driving and of the WBL, schemes based on the tra- 
ditional approach Q, but working directly in the time 
domain have been put forward [H, [l^- AH time con- 
volutions, which result from the projection onto the de- 
vice states, are transformed into matrix-matrix multipli- 
cations using a time-discretization scheme. In contrast, 
the formalism presented in Refs. [3, [3 is based on prop- 
agating the wave function of the full system (device and 
reservoirs). This is accomplished by using the Cayley 
propagator and then projecting on the subsystem of in- 
terest. This formalism provides a natural way to work 
within the so-called partition-free approach, where the 
time-dependent voltage pulse is considered to act on the 
total system. 



In this article we present a general propagation scheme 
which is also based on the NEGF formalism Q and al- 
lows to obtain device-related observables and the electric 
current as a function of time. Our formulation, how- 
ever, relies on a set of coupled equations of motion for 
quantities with only one time argument. In this way we 
can avoid time convolutions and standard methods for 
integrating the equations may be applied. As with the 
standard formulation, the key issue consists in perform- 
ing integrals over reservoir states which eventually lead 
to tunneling self-energies. In this context we propose us- 
ing an auxiliary-mode expansion which allows to treat 
finite temperatures. We provide a numerical implemen- 
tation of our scheme for two relevant cases — the wide- 
band limit and a level-width function given by a sum 
of Lorentzians. The methods are applied to transport 
through a randomly fluctuating level and the transient 
response of a DQD to a voltage pulse. 

In the remaining part of the introduction we briefly 
discuss the general setup [Sec. II A| and then repeat the 
findings of the standard NEGF formalism in the context 
of our propagation scheme [Sec. II Bj . 

A. Setup 

We take the usual threefold setup consisting of a device 
(system) which is coupled to two electron reservoirs. The 
coupling is due to tunneling through a barrier. The total 
Hamiltonian is 

H = Hu + Hn + i/DR . (1) 

The device is described in terms of discrete energy levels 
e„(t) which may be coupled through Vnmit), 

Hu=Yj £n{t)cicn + Vnm{t)clc„, . (2) 

n n^m 

The operators {cj^} and {c„} denote the creation and 
annihilation of an electron in state n. The reservoirs are 
described by non-interacting electrons and the Hamilto- 
nian reads 
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with single-particle energies of the form eak{t) = £°j, + 
^akit)- Finally, the coupling Hamiltonian is 

= E E (^) ^ikCn + h.C. , (4) 
n ak 

with {T^nl denoting the couplings between device and 
reservoir a — L, R; {6^^.} and {bak} are electron creation 
and annihilation operators for reservoir states, respec- 
tively. 

Regarding the time dependence of the reservoirs and 
the device we adopt the Caroli partition scheme p^ . 
i.e. all sub-systems are separated at t = —oo and in 
their respective equilibrium state. Any time dependence 
only sets in after eventually coupling the different parts. 
Consequently, the single-particle occupation probability 
in the reservoirs is determined by e^;,; the time depen- 
dence Aak (t) of the reservoir energies appears as a phase- 
factor only. The situation where the chemical poten- 
tials and therefore the occupation probabilities are time- 
dependent has been critically discussed before Q- 



which applies to Green functions as well as to self ener- 
gies, one can find advanced (retarded) self-energies in Eq. 
([5]) in terms of greater and lesser functions. The level- 
width function Fq in Eqs. ^ depends on the density of 
states Pa{£) of reservoir a and the coupling Ta^n{£) of 
device level n and the reservoir state at energy e, 

[ra(e, ti, t)]„„ =2^Pa(£)T„,„(£, t)T:^^{e, h) 

t 

X exp{i J dt2A^{e,t2)} . (8) 
ti 

Replacing the advanced and retarded quantities in Eq. 
([S]) by using Eq. ([7]) one can rewrite the expression for 
the current in a very compact form, 

Jc{t) = 2eReTr {n„(t)} . (9) 

The current matrices Tla{t) are given by the following 
expression 



B. Time-Dependent Current and Non-Equilibrium 
Green Functions 

By applying the Keldysh formalism to non-equilibrium 
Green functions it is possible to obtain a general formula 
for the time-dependent current in the setup introduced in 
Sec. II Al The current Ja through the barrier connecting 
lead a and the device is given by 0, Q 

t oo 

Jait) = 2eReTr i [ d^i [G<(t, t) 



+G'-(t,ti)S<(ii,t)] I . (5) 

Throughout the paper we adopt units with h ~ 1. In Eq. 
([5]) G< and G'' are lesser and retarded Green functions 
and and are advanced and lesser self-energies, re- 
spectively. All boldface quantities are matrices related 
to the device states, e.g., G<(t,ti) = G^^{t,ti). Prod- 
ucts are to be understood as matrix multiplications. The 
greater and lesser self-energies are explicitly given by 

S>(ii,i) ^ -ij ^/„(e)e--(*i-*)r4£,ti,i)(^a) 

S<(ti,i) = i J ^/„(e)e--(*-*)r„(£,ii,t).(6b) 

As indicated at the end of the previous section the Fermi 
distribution, fais) = /(/5(e— /^q)), characterizes the equi- 
librium state of reservoir a with the chemical potential 
/ia and inverse temperature f3 — (k-oT)^^ at io ~ — oo. 
It is /^(e) = 1 — /a(e)- Using a relation for two-time 
functions, 

x--(t,t') = ±e(±tTO {x>{t,t!)-x<{t,t!)\ , (7) 



n„W = J dt2{G>{t,t2)-S<{t2,t)-G<{t,t2)i:>{t2,t)) , 
to 

(10) 

where the first and the second term describe electrons 
tunneling into and out of the device, respectively. Equa- 
tions (O and pO)) have been discussed in the context 
of current conserving self-energies Q. The conceptually 
new approach presented in this article consists in con- 
sidering Tla{t) as an independent entity. In particular, 
opposed to correlation functions such as G^{t,t2) the 
current matrices Tla{t) only depend on a single time ar- 
gument. Therefore, they are fully determined by a single 
equation of motion. This circumstance provides the basis 
of our propagation scheme, which is presented in Sec. [TTl 
Moreover, in order to calculate the expectation value of 
any device observable Od it is advantageous to use the re- 
duced single-electron density matrix, a{t) = Im G< {t, t). 
The expectation value is then given by 

(OD(i)) -TrD{ODT(t)} . (11) 

Similar to the current matrices 11^ {t) the density matrix 
only depends on a single time argument and one has the 
following equation of motion 

i|^(t) = mt),rT{t)]_ 

+i^(n„(t) + ni(i)) , (12) 

which depends on the current matrices !!„ (t) . The bold- 
face Hamiltonian H(t) = [Hd]^^ is obtained from the 
device Hamiltonian in Eq. ^ . Equation (fT2|) is found by 
using G'^{t',t) — — [G^(i,t')]^ and from the equations 
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of motion for greater and lesser Green functions and 
G<, 

i|-G^(t,t') = H(i)G^(<,0 



dt 



J dt2'Ei,it,h)G^it2,t') . (13) 



The total self-energies are sums of the tunneling 

self-energies for each reservoir. The Green functions may 
also be obtained from the Dyson series leading to an in- 
tegral equation Q . 



II. CURRENT MATRICES AND AUXILIARY 
MODE EXPANSION 

In order to arrive at a viable propagation scheme we 
will rewrite the equations of motion given above by intro- 
ducing energy-resolved quantities. This form allows for 
applying an auxiliary-mode expansion which replaces the 
energy integrals by finite sums. The number of (matrix) 
equations to be propagated is determined by the size of 
the expansion. 



A. Energy-Resolved Current Matrices 

First we assume factorizing momentum and time 
dependence of the tunnel coupling, Ta^ni£,t) — 
Ta^n{£)ua,n{t)- The Same ansatz has been used in Ref. 
Q for the non-interacting resonant-level model. For no- 
tational convenience we consider in the following only the 
case of a common time dependence of the coupling for all 
device states, i.e. Ua^n{i) = Ua{t). Equation ([H]) becomes 



Ta{e,ti,t) = M*(ii)ua(i)rQ(£)exp{i J dt2Aa{e,t2)} ■ 

ti 

(14) 

Next, we define energy-resolved self-energies as 

S>(e;ti,i) - -i<(ti)/„(e)e--(*^-*)r„(£) 

t 

xexp{i J dt2Aa{e,t2)} , (15a) 
ti 

S<(e;ti,i) = i<(ti)/„(£)e-'^(*i-*)r„(e) 

t 

xexp{i J dt2Aa.{eM)} ■ (15b) 



In terms of these expressions the full self-energies are 
given by 

S|(ti,0 = u^{t) I d£S|(£;ii,i) , (16) 

which follows from Eq. (|6|). Using the definitions above 
we introduce energy-resolved current matrices, 

t 

n„(e;t) = J dt2 {G>{t,t2)-S<{e;t2,t) 



-G<{t,t2)-S>{e;t2,t)) 



(17) 



From Eq. ([TT]) one finds Ila{e;to) = 0. The expression 
for the current given by Eq. ^ becomes 

Ja{t) = 2eReJ2uc.it) J den„,„„(e;t) . (18) 



Therefore, the diagonal elements Tla.nni^'it) may be in- 
terpreted as the current flowing from the reservoir state 
at energy e to the system state n. The total current 
through the barrier is then given by the sum of all pos- 
sible currents. 

The equation of motion [Eq. p2| ] for the reduced 
single-electron density matrix cr of the device becomes 



i|^(i) = [H(<),^(i)]^ 



(19) 



J de(u^{t)U^is-t)+ul{t)Ui{e;t) 



which now contains the energy-resolved current matrices. 

Due to the definitions of the energy-resolved self- 
energies, their time derivatives, 



d_ 

dt 



S|(£;ti,i) =i(e + A,(£,i))S|(e;ii,i) , (20) 



and by using Eq. (jl3p . one gets an equation of motion 
for the energy-resolved current matrices, 

i§-U^ie;t) = -^<(f)(^(t) -/„(£)) r„(e) 

Ot ZTT 



-{H(<)-(£ + A„(e,t))}n„(e;t) 
^<,(0 J de'Vlo.a.'{e,e'-t) , 

Oi' 

(21) 



where a new quantity ^aa' has to be introduced. It 
contains all contributions from the time derivative of the 
greater and lesser Green functions, which give rise to a 
double time integral. Consequently, its definition is 
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n„„,(e,e';i) = J dt2 J dhi:l,ie'-,t,ti) [G>ih,t2)'E<{e;t2,t) - G<{h,t2)'E>{e;t2,t)] (22) 
J dt2 J dil [S<(e';i,ii)G"(<i,t2)S>(£;i2,i)-S>(e';i,ti)G"(<i,t2)S<(£;t2,i)] • 



^0 io 



We replace the retarded self-energies and the advanced Green function again using Eq. ([7]), but instead of showing 
the result we rather give the equation of motion, which is easily obtained from Eq. p2p , 



. d 



+ {(e' + A„,(£',i))~ (e + A„(e,t))} (£,£'; i) , 

I 



(23) 



with the initial conditions fiacf (e, e'; to) = 0. The equa- 
tions of motion given by Eqs. (fT9|l . ([2T|l and ([23]) provide 
a closed description of the non-equilibrium dynamics of 
the device. A similar set of equations has been found 
recently [iTj . where it was derived from a hierarchy for 
the many-body density matrix. The identification of 



-iHn and 



(24) 



renders their equations identical to the ones given 
above. This provides an independent verification of the 
density-matrix approach [13] and shows that the hier- 
archy derived therein yields the exact dynamics for non- 
interacting electrons under the assumptions stated above. 

The full single-particle density matrix has a size of 
(A^D + ^r)^, where N-£, and A^r are the number of 
single-particle states in the device and the reservoirs, 
respectively. In the present case we have to propagate 
X (iVfi + 1)^ quantities with 11 and 11^ counting inde- 
pendently. Therefore, the complexity of Eqs. p^ . (PT|) 
and is at least the same compared to calculating 
the full single-particle density matrix. In particular one 
has to deal with a continuum of states and consequently 
the utility of the method depends on finding an efficient 
strategy for performing the energy integral. In the fol- 
lowing subsection we will provide such a method based 
on the expansion of the Fermi function and making use 
of the residue theorem. The same strategy has been suc- 
cessfully applied to the propagation of non-Markovian 
quantum master equations involving bosonic [Tsj and 
fermionic reservoirs H^, fioj . The formulation in terms 
of energy-resolved quantities depending on a single time 
argument turns out to be beneficial in this context. In 
order to propagate each matrix only the value of the pre- 
vious time step has to be known. References to past 
times [l^, [13] are not necessary. This comes at the cost 
of having to propagate the two-energy quantity flaa'- 
However, as we will show in the next section one can ef- 
fectively reduce the associated numerical costs by using 
an auxiliary-mode expansion. 



B. Auxiliary-Mode Expansion 

The general idea of the auxiliary-mode expansion con- 
sists in making use of contour integration and the residue 
theorem. To this end the Fermi function is expanded in 
a sum over Np simple poles. 



(25) 



with Xap = l^a^Xpl (3 and IvaXp > 0. The well-known 
Matsubara expansion 20] is an example for such a de- 
composition. Its major disadvantage consists in a poor 
convergence behavior especially for low temperatures. A 
particular efficient alternative is presented in appendix 

El 



1. Wide-Band Limit 

As a first application we consider the WBL, i.e. 
rQ(e) — const. From the definition of the self-energies 
(O and the expansion of the Fermi function [Eq. (|25p ] 
one obtains for t > ti, 



(26) 



+«„(t)^ir„<(ii)e'^i'^*^^-(*^), 
p 

where xtp{t) = Xap + ^a{t)- Analogously, one finds for 
the lesser self-energy 



(27) 



Thus, the expansion of the Fermi function leads to an 
expansion of the self-energies into a sum of exponentials. 
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Due to the WBL one also gets one term proportional to a 
delta function. We introduce auxiliary self-energies Sap, 
which incorporate the exponentials, i.e. 



J-a{t) '^ap{tl, t) , 



(28a) 



= ir„<(ti)e'A'**«-(*^' , (28b) 

which implies S„p(i,t+) = jjTau'^{t). Next, we insert 
the expanded self-energies into the definition of the cur- 
rent matrices (ITOl). 



(29) 



Tla{t) = i|u„(t)P(l-2cr(t))r„ 



and obtain an expansion in terms of auxiliary current 
matrices, 

t 

n„p(t) = J dt2{G<{t,t2)^apit2,t) 
to 

-G>(t,t2)S„p(<2,i)) . (30) 
Their equation of motion is easily found. 



(31) 



H(<) - -r{t) - xUt)i ] u^p{t) , 



where T{t) = J2a' Wa' {tjl'^'i' a' ■ The coupled equations 
of motion ([T2|) and ([3T|) allow with Eq. ((29|) for a com- 
plete description of the non-equilibrium dynamics of the 
device. Comparing Eqs. (pij) and ([?T|) suggests that 
S^ap,Q'p'(i) = — 5Ua'(0r'a'nQp(t)(5pp'. Tlius, an addi- 
tional equation of motion for is not needed for the 
WBL. 



2. Lorentzian Level-Width Function 

The next application we consider is the case of a 
Lorentzian level- width function (LLWF). We take a gen- 
eral ansatz of the form 



ra(e) = E 



■ ai 



■ Oil 



£-£al- iWae iWat 



(32) 



with Wai > and T^^ = T^^aeWae- Equation ^ 
might be used as a par ametrization of an arbitrary level- 
width function p^. Il9||. Now, we can plug Eq. into 
the definition of the self-energies [Eq. ^] and evaluate 



the energy integral by means of contour integration. This 
procedure yields for t > ti 



"P„g-i(e„f+iM^af)(ti-t) 



E^r„(xip)e-' 



X exp{i J dt2Aa{e,t2)} 
ti 

-E^r„(xip)e---(*-*) 



(33a) 



t 

exp{i J dt2Aaie,t2)} 



(33b) 



where /^^ = fa{£ai + iW^^) indicates that the expansion 
given in Eq. ((25)) should be used to calculate the Fermi 
function at the position of the pole i. The self-energies 
are thus given by a finite sum with Ni^ -\- Np terms. For 
convenience we combine the two indices p and i yielding 
a single index x = {£,p}. The coefficients and exponents 
are combined in a similar way, 

{±r± /„(£„, ±iTy„,),±ir„(x±p)}(34a) 

xtx = {eoii±WM,xtp} ■ (34c) 

Using these conventions the self-energies can be written 
in a compact form, assuming t > ti we have 



S|(ii,t) = u„(t)Es|,(ti,t) 



(35a) 



Si.(ii,<) - <(ti)rii+e'A'^*«i^(*^) , (35b) 

where xt^xi^) — xtx + Aa{t)- The auxiliary self-energies 
S<^ are simply exponentials. The respective auxiliary 
current matrices can be calculated in analogy to the 
energy-resolved current matrices, i.e. 

t 

na,(t) = J dt2{G>{t,t2)-E<^{t2,t) 
to 

-G<it,t2)i:>At2,t)) . (36) 
Their equation of motion is then given by 

i^n^xit) = ui{t)T<^+ + ui{t)cT{t) (r>i+ - r<i+) 
+ {ii[t)-xtx{t)) ^Ut) 
+ E<'W" 

(t) . (37) 
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The initial condition Tlaxito) = follows from Eq. ([36)1 . 
Notice the similarity to the energy-resolved current ma- 
trices given by Eq. ([2T|) . In particular, we also have a two- 
mode quantity flax.a'x' appearing in the equation of mo- 
tion. Its definition is again in full analogy to the energy- 
resolved case given in Eq. (HID, but with S<(e';t,ii) 

replaced by S^,^,(t, ii). Also the equation of motion is 
similar to the energy-resolved case [Eq. ([23]) ] , 



i^a'W (r>g-r<;-)n„,(t)(38) 
+inl^,{t) (r>i+-r<i+)<(t) 
+ {xa'x'it)-xtx{t)) n^.,a'x'{t) . 

At < = io one finds ^ax,a'x'{to) = 0. It is interesting to 
notice that for x — p,x' = p', i.e. both indices represent 
an auxiliary mode resulting from the Fermi-function ex- 
pansion [Eq. (|25|l]. one gets ^i^ap^a'p'it) OC 0,ap,a'p'{t). 

Taking the initial condition into account it follows that 
^ap.a'p'{i) = for all times. This is consistent with 
the energy-resolved expression [Eq. ([121)] where any ref- 
erence to the Fermi function is absent. Consequently, in- 
stead of propagating (A^l + A^f) x (A^l + A^f) J^-matrices 
we only need to consider A^l x (A^l + '^Np) matrices for 
each reservoir index a. Since typically A'l ^ A'p the 
memory requirement of the proposed method scales with 
A^L X A'f and the computational time requirement scales 
with A^T X A^L X A'f , where A^t is the number of time steps. 
Notice that in spite of having to use two-energy quanti- 
ties, using the auxiliary-mode expansion for the Fermi 
function yields a scheme, which scales linearly with the 
number of modes and thus allows for a particularly effi- 
cient propagation. 



III. APPLICATIONS 

We apply the proposed propagation scheme to two sit- 
uations: a resonant-level model with a randomly fluctu- 
ating energy level and a DQD system driven by finite 
bias-voltage pulses. These two situations demonstrate 
that our scheme is especially suited to study a strongly 
fluctuating driving and realistic experimental pulses in- 
cluding structured reservoirs. 



A. Fluctuating Energy Level 

As a first application we consider a resonant-level 
model with a single randomly fluctuating energy-level 
SdU), which is given by a Gaussian stochastic process 
[28|. An analytic expression for the current is given in 
appendix |B] 

The device Hamiltonian [Eq. ^] is simply. 



(39) 



with "d" denoting the single-electron state of the device. 
All matrices become scalars and the respective equations 



of motion are scalar equations. The stochastic process 
Edit) is fully characterized by the first and second mo- 
ments. 



(edit)) = 0, 
(ed(i)ed(i')> = c{t-t') 



(40a) 
(40b) 



Here, we take £d(i) as realization of an Ornstcin- 
Uhlenbeck (OU) process, which yields for the correlation 

2 

function c{t — t') = ^ exp[— «(< — *')]. The OU process is 
characterized by two parameters, the inverse correlation- 
time K, and the noise amplitude rj f2l'|. 

Considering the WBL and using a symmetric coupling 
to the left and right reservoir, Fl = Fr = F/2, we sud- 
denly connect the device and the reservoirs at i = 0. 
The reservoirs are further characterized by chemical po- 
tentials = 2F, /xr — F and temperature k^T — O.IF. 
Thus, without stochastic driving the energy level is not 
located in the transport window and a non- vanishing cur- 
rent is a result of the broadening due to the coupling to 
the reservoirs. 

The equations of motion obtained in Sect. Ill Bl are 
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FIG. 1: Time-resolved occupation (A) and net-current (Jnet) 
for different values of the noise amplitude. Noise averages are 
obtained from 20000 realizations and n = O.SF. The arrows 
indicate the time-averaged values ( Jnct) obtained by sampling 
the current for times t > lO/F. 
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propagated using a weak second-order Runge-Kutta 
scheme [l^l with a constant time step ^] 6t = 0.01/r. 
We use A'f = 240 auxihary modes for all calculations. 
The resulting time-resolved occupation, N(t), and net 
current, Jnct{t) — [Jhit) — JB.{t)] /2, are averaged over 
20000 realizations of the stochastic process. Figure [T] 
shows the averages {N{t)) and {Jnct{t)) for k = O.SF and 
three selected values of i] = 0.5, 1.0, 3. OF. We also show 
the case of no stochastic driving. One sees a transient 
response to the sudden coupHng for times t = . . . 10/F 
and the eventual settling to a stationary value. In all 
cases shown in Fig.[T]the stationary current is larger than 
for the case without any noise; but the dependence on rj 
is non-monotonic. 

In order to quantify the stationary current we take 
the time average (Jnot) for the time interval starting at 
t — lO/F. Figure shows the obtained time-averaged 
current as a function of the noise strength rj and for 
various values of the inverse correlation-time k. The 
time-averaged current exhibits a pronounced maximum 
as a function of noise strength; the transport through 
the energy-level is stochastically enhanced. This ef- 
fect reminds of the phenomenon of stochastic resonance 
[2^. The observed behavior is a result of additional 
broadening due to the stochastic driving 'S|]. The cur- 
rent is proportional to the area under the spectral den- 
sity, A{e) = —2lmG'^{e), within the transport window 
given by {/iL,/iR}, cf. appendix IB] For increasing noise 
strength the spectral density becomes broader and has 
more weight in the transport window. However, the 
height of A{e) decreases at the same time which even- 
tually leads to a decrease in the area in the transport 
window. These findings are corroborated by the analyti- 
cal result [Eq. (jB4|) ]. which is also shown in Fig. [21 The 
numerical results agree very well with those results. 
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FIG. 2: Time-averaged current vs noise amplitude for dif- 
ferent values of the inverse noise correlation-time k. Noise 
averages are obtained from 20000 realizations. Error bars in- 
dicate 95% confidence interval for the sample mean. Full lines 
denote the analytical result given by Eq. (|B4|l . 



B. Double Quantum Dot 

As a second application we will now discuss the re- 
sponse of a DQD to a voltage pulse. The device consists 
of two QDs which are coupled in series. Each dot is also 
coupled to an electron reservoir. This setup resembles a 
typical experimental situation (see for example [l|). 

The DQD is modeled by a two- level system, i.e., one 
localized energy- level per dot. The device Hamiltonian 
[Eq. (121)] is then. 



Hr 



£d{t) c\cd + Vc\c^ + h.c. 



(41) 



with "1" and "r" the localized single-electron states. The 
time-dependent bias- voltage is assumed to act on the en- 
ergies in the following way: AL(t) = — AR(i) = Vbias(i)/2 
and £\{t) = —£r{t) = Vbias(i)/4. Initially, the chemical 
potentials /^l and /iR and the QD energies ei^r are zero. 
The temperature is k^T — O.IF for both reservoirs. Since 
the two dots are coupled in series, the level-width func- 
tions contain one non-zero element, 



F(e)/2 






F(£)/2 



The matrix element F(e) is either constant in the case of 
WBL, or is taken to be a single Lorentzian [23 |. 



F(e)=F- 



1^2 



(42) 



The latter is compatible with the general ansatz given 
in Eq. ((32|) and is chosen such that WBL is attained for 
IF — > 00. For the time dependence of the bias voltage we 
take a rectangular pulse, i.e.. 



Vbias(i) 



tanh ( — 



tanh 



, (43) 



which is characterized by the pulse length t-p. The finite 
switching time ts reflects the experimental situation (e.g. 
as reported in Refs. [H, Q). In the following calculations 
we use is = 1/F and Vmax = 3F. The equations of motion 
obtained in Sect. Hi] for the WBL and the LLWF, respec- 
tively, a re p ropagated using a fourth-order Runge-Kutta 
scheme |25| with constant time step 5t — 0.02/F. We use 
N-p = 120 auxiliary modes for all calculations. 

Figure [3] shows the numerically obtained current Jl as 
a function of time t for different widths W in response 
to the same pulse of length t-p = 20/F. The current 
shows a transient behavior at the beginning and after 
the end of the pulse. For sufficiently long pulses it set- 
tles to a new stationary value according to the plateau 
bias voltage Vbias — Knax- Notice that this situation for 
a structured reservoir is different from initially having 
fJ-h ^ fJ-R = Knax- In the latter case the chemical poten- 
tial and the center of the level-width function [Eq. (|42p ] 
are shifted with respect to each other. The two distinct 
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FIG. 3: (color online) Time-resolved current through left bar- 
rier Jl for different values of W in Eq. (|42p driven by a bias 
voltage pulse according to Eq. (|43p with Vbias = SF, ts = l/F 
and tp = 20/r. The WBL corresponds to ^ oo. 



situations are illustrated in Fig.lH We adopt the physical 
relevant situation shown in the right panel (see also Ref. 

0). 

At any rate, the stationary current is found to van- 
ish for T4^ ^ 0, which is an artifact of the level-width 
function given by Eq. (|42|) . The ringing behavior at the 
beginning and after the pulse is qualitatively similar for 
all values of W. However, for small W the damping of the 
current oscillations is weaker. In experiments the direct 
observation of the ringing may be obscured by capacitive 
effects or the resolution of the ampere meter. Therefore, 
one considers the time-averaged or time-integrated cur- 
rent as a function of pulse length [1, [2^ . The latter yields 
the number of pulse-induced tunneling electrons. 



Np{tp) = / dt [J{t) - Jo] 



(44) 



where Jq is the stationary current without pulse and 
J{t) = Jh{t). In Fig. [5] we show A^p as a function of 
pulse length tp for various values of W. One observes 
an increase in the number of tunneling electrons with 



3 



r 























n 






r 



FIG. 4; (color online) Energy scheme for W = F at the 
plateau of the pulse with Vbias ~ Vmax ~ 3F. Outermost 
parts show /Q(e)r(e) for a = L, R (blue/dark-shaded areas). 
The inner part shows spectral densities of left and right levels. 
Left panel: Vbias = Ml ~ MR ^-nd Al — Ar — 0. Right panel: 
Vbias = Al - Ar and ml = MR = 0. 




FIG. 5: Number of pulse-induced tunneling-electrons A'^p vs 
pulse length tp. Symbols indicate numerical results for dif- 
ferent values of W in Eq. (I42|l . Straight lines show result of 
linear fit in the respective range. The dashed line gives the 



two intercepts Np 
WBL. 



and tp, respectively, for the case of the 



increasing pulse length. Remembering the time depen- 
dence as shown in Fig. [3] it is clear that for short pulses 
A'p is dominantly determined by the transient part of the 
current. For sufhciently long pulses, however, the main 
contribution comes from the new stationary current Jstat 
and one expects A^p cx tp with a slope given by Jstat • This 
asymptotic behavior is shown in Fig. [5] by the straight 
lines which have been obtained from a linear fit to the 
numerical data. The slope was fixed by independently 
calculating Jstat using stationary NEGF formalism @]. 
The fitting procedure yields the A'p-intercept denoted by 
Np (cf. dashed line in Fig. O, which provides a mea- 
sure of how transient the current response actually was. 
If the current would instantaneously switch to the new 
stationary value one would get A^p = Jstat^p and the A^p- 
intercept would vanish. Non-vanishing values of A^* re- 
flect the additional transient contributions to the current. 
Figure [5^ shows a stronger transient response for smaller 
values of W which is in accordance with the observations 
for the time-resolved current. In any case the net ex- 
cess is positive since the transient response following the 
switching-on outbalances the one after the switching-off. 

Using the A^p-intercept and the stationary current one 
can also calculate the pulse length that would be neces- 
sary to yield the same number of tunneling electrons if 
the DQD would switch instantaneously, t* = N*/ Jstat 
(cf. dashed line in Fig. [5]). This quantity is shown in 
Fig. [61d. It gives a measure for the pulse length at which 
transient and stationary contributions are of similar size. 
Therefore, the transient response for pulses with tp 3> t* 
becomes negligible. 
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IV. SUMMARY 

We have presented a propagation scheme for time- 
dependent electron transport which is based on non- 
equilibrium Green functions. It relies on quantities with 
a single time argument which allows for a straightfor- 
ward numerical implementation with standard differen- 
tial equation solvers. 

The basis of our scheme is a reformulation of the well- 
known expression [3| for the current J{t) by means of 
the density matrix (T{t) and newly introduced current 
matrices n(t), cf. Eq. (fTU|) . Decomposing these matri- 
ces into energy-resolved expressions allows to obtain a 
closed set of coupled equations of motion for (T(t) and 
Il{e,t). Thereby, one has to consider another energy- 
resolved quantity fl{e,e';t) given in Eq. ([^ . The equa- 
tions of motion are given by Eqs. p9)) . ((2T|) and ((23|) . 

For a numerical implementation of these equations we 
propose using an expansion of the Fermi function and 
a parameterization of the level-width function by a set 
of Lorentzians (Tsj . The error made by truncating the 
expansion can be reduced by applying a fast converging 
decomposition (27j . The matrix equations to be solved 
are p2|) , p7|) and p8|) , respectively. In the often applied 
wide-band limit the set of equations simplifies since $7 
can be found explicitly in terms of the current matrices, 
cf. Eq. (EID. 

Finally, we have applied our scheme to two illustrative 
examples: the randomly fluctuating energy level and the 
response of a DQD to a voltage pulse. In both cases 
a non-trivial driving was involved. For the DQD we 
demonstrated the influence of structured reservoirs on 
the transient current response. This transient contribu- 
tion may be quantified by using the number of pulse- 
induced tunneling electrons [Eq. p4)l ] as a function of the 
pulse-length. For the fiuctuating energy level we showed 
a good agreement of our numerical calculations with an- 
alytic results obtained for the stationary current. More- 
over, we found an enhancement of the current due to 
the stochastic driving. The study of this effect in more 
complex systems might lead to interesting new applica- 
tions. In general, we expect our method to be a valuable 



W [T] 
4 2 1 




W [T] 
4 2 1 




tool for investigating time-resolved electron transport in 
nanoscale devices. 
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APPENDIX A: EXPANSIONS OF 
SELF-ENERGIES 

In order to perform the energy integration in Eqs. 
we expand the Fermi function in terms of a finite 
sum over simple poles. This procedure yields the ex- 
pression given in Eq. (|25p . The poles are given by 
xtp = tJ'a+Xp/6 = (xapY ■ Instead of using the Matsub- 
ara expansion [20|, with poles Xp = i7r(2p— 1), we use a 
partial fraction decomposition of the Fermi function [l^ , 
which converges much faster than the standard Matsub- 
ara expansion. Furthermore, it allows to estimate the er- 
ror made by truncating the sum [Eq. 1^5]) ] at A^f terms. 



For this decomposition the poles Xp = ±2^/z^ are given 
by the eigenvalues Zp of the N-p x N-p matrix [23| 



(Al) 



We take the branch of the root y^z^ such that Im (xp) > 
for all p. Thus all poles Xp iXp ) ^-re in the upper (lower) 
complex plane. 

Given the expansion [Eq. (PS)) ] one can evaluate the 
energy integrals by a contour integration in the upper 
or lower complex plane depending on the sign of t — ti . 
Thereby, the integration becomes a (finite) sum of the 
residues. 



APPENDIX B: NOISE-AVERAGED CURRENT 
FOR RLM 

The noise-averaged net-current, (Jnct(i)) = 
{Jh{t) — Jnit)) /2, can be obtained from the gen- 
eral expression for the time-dependent current [Eq. 



(Jnct(t)) 



e Re Tr 



[S<(ti,<)-S<(ii,t)]} , (Bl) 



0.5 1 1.5 
1/W [ l/T] 



where a symmetric coupling, Ti^{e,ti,t) = T^{e,ti,t), is 
assumed. For the resonant level model all quantities are 
scalars and in particular for the setting considered in Sec. 
nil Al one has 



FIG. 6: Results of the linear fit to Np{tp) shown in Fig. (5] a) 
A^p-intercept and b) fp-intercept. 



G''{t,ti) = -ie(t-ii)exp 
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(B2) 

In order to evaluate Eq. (|B1[) we need the average of 
the fluctuating exponential function in G^{t,ti) which is 
obtained by using the cumulant expansion, i.e., 



exp 



exp 



exp 



exp 



dt'£d{t') 



dri / dT2 (ed(Ti)£d(T2)) 



t 2- 



c(w) 



div 
duj c{uj) 



ti 
4 sin 



2 / - ti) 



.(B3) 



In the derivation we have used the properties of the noise 
[Eqs. (501)] and introduced the Fourier transform of c(ti — 
T2) which is denoted by c(w). 

Thus, the noise-averaged retarded Green function does 



only depend on the time difference and the time-averaged 
current is given by a Landauer-type expression Q 

where the spectral density A{e) is given by the time and 
noise averaged retarded Green function. 



A{e) = -2Im j dr {G'{t,t - t)) e'"^ 


^^gier-r|r|/2 



(B5) 



exp 



i/^^4sin^f^ 
2 / 27r w2 \ 2 



For white noise one has c(w) = 7 and the fluctuations 
lead to a trivial broadening of the spectral density. 
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